Point defects in two-dimensional colloidal crystals: simulation vs. 

elasticity theory. 
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Using numerical and analytical calculations we study 
the structure of vacancies and interstitials in two- 
dimensional colloidal crystals. In particular, we com- 
pare the displacement fields of the defect obtained nu- 
merically with the predictions of continuum elasticity 
theory for a simple defect model. In such a compari- 
son it is of crucial importance to employ corresponding 
boundary conditions both in the particle and in the con- 
tinuum calculations. Here, we formulate the continuum 
problem in a way that makes it analogous to the elec- 
trostatics problem of finding the potential of a point 
charge in periodic boundary conditions. The continuum 
calculations can then be carried out using the technique 
of Ewald summation. For interstitials, the displacement 
fields predicted by elasticity theory are accurate at large 
distances, but large deviations occur near the defect for 
distances of up to 10 lattice spacings. For vacancies, 
the elasticity theory predictions obtained for the simple 
model do not reproduce the numerical results even far 
away from the defect. 

1 Introduction 

Many properties of crystalline materials are strongly af- 
fected by the presence of imperfections in the crystal 
lattice. In particular, point defects such as vacancies 
and self interstitials have a profound influence on the 
mechanical, optical, and electrical behavior of the ma- 
terial. Recent advances in experimental techniques for 
the manipulation and observation of colloidal systems 
[H [2] now permit to study the fundamental properties of 
point defects in condensed matter systems with "atom- 
istic" space and time resolution. Using optical tweezers 
to manipulate individual colloidal particles, Pcrtsinidis 
and Ling [21 [H [5] have generated point defects in two- 
dimensional crystals and have studied their stable struc- 
tures, interactions and diffusion. In other experimental 
work, Maret, Griinberg and collaborators [SI El [5] have 



investigated the effective interactions of thermally ex- 
cited topological defects in crystals of paramagnetic col- 
loidal particles and discussed the significance of these 
interactions for 2d-melting, which according to the cele- 
brated Kosterlitz-Thouless-Halperin-Nelson- Young the- 
ory [9] , involves the formation and dissociation of topo- 
logical defect pairs. Point defects also play an important 
role in the two-dimensional electron lattice, the so called 
Wigner crystal [TU], in which they carry implication for 
the melting mechanism [TTJ [T^ [T3] , and for the the con- 
jectured supersolid phase of Helium 4 [M], in which case 
the attractive interactions of vacancies and interstitial 
may lead to expulsion of defects from the crystal thus 
preventing formation of a supersolid [15j . 

From experiments O [H [5] and computer simulations 
|16l 117] it is known that vacancies and interstitials in 2d 
colloidal crystals can occur in various stable configura- 
tions with symmetries that differ from the symmetry of 
the underlying lattice. In the present article we study 
the structure and energetics of such point defects in a 2d 
crystal of soft spheres using computer simulations and 
analytical calculations. In particular, we address the 
question of how accurately the disturbances created by 
point defects can be rationalized in terms of elastic con- 
tinuum theory. Due to the long range nature of elastic 
displacement fields, in carrying out such a comparison 
it is critical to use corresponding boundary conditions 
in the particle and continuum calculation. Similar pe- 
riodic image effects due to elastic interactions need to 
be taken into account also in the atomistic modeling of 
dislocations [HI [19] . As we show below, the structure of 
point defects in a system with periodic boundary condi- 
tions can be determined within elasticity theory with the 
technique of Ewald summation familiar from the com- 
puter simulation of systems with electrostatic interac- 
tions [20] . This technique has been used before to adapt 
the interaction of dislocations to periodic boundary con- 
ditions [ini [211 [22 • Here we use Ewald summation to 
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solve the equilibrium condition of elasticity theory and 
calculate the displacement field of a simple point defect 
model under periodic boundary conditions. While elas- 
ticity theory accurately describes the lattice distortion 
caused by point defects in the far field, non-linearities 
and discrete lattice effects dominate the defect structure 
near the defect. 

Although all the numerical studies discussed in this 
paper are carried out for two-dimensional crystals of soft 
spheres, simulations performed for three-dimensional 
crystals of various structures and with different interac- 
tion potentials, including Gaussian core, Lennard- Jones, 
screened electrostatic, and l/r'^-interactions, indicate 
that the phenomena described here are common to many 
atomic and colloidal systems. 

The remainder of the paper is organized as follows. 
In Sec. [2] we describe how wc determine the displace- 
ment fields of point defects numerically, discuss whether 
such calculations should be done at constant pressure or 
at constant volume, and present the displacement fields 
caused by interstitials and vacancies in various configu- 
rations. The elasticity theory formalism we use to an- 
alyze the displacement patterns of point defects is de- 
veloped in Sec. [3l In this section, we also discuss the 
analogy between elasticity theory and electrostatics that 
enables us to use the method of Ewald summation to 
obtain displacement fields from elasticity theory. These 
displacement fields are compared to those obtained nu- 
merically in Sec. [51 A summary and conclusions are 
provided in Sec. [5l 

2 Displacement fields 

Throughout this paper, we use the Gaussian core model 
as a generic model for a system of soft spheres [26l [27l 
I28j . In this purely repulsive system, pairs of particles 
interact via 

z)(r) = eexp(-r^/cr^) (1) 

where r is the inter-particle distance and e and a set 
the energy and length scales, respectively. In the follow- 
ing, energies are measured in units of e and distances in 
units of (J. The Gaussian core model, often studied in 
soft condensed matter science, accurately describes the 
short-ranged effective interactions between polymer coils 
in solution [29| . Depending on temperature and density, 
the three-dimensional Gaussian core model can exist as 
a fluid, a bee- or an fee-solid [27]. In two dimensions, the 
perfect triangular lattice is the lowest energy structure 
at all densities [30]. Since Gaussian core particles are 
purely repulsive, they can form stable crystals only at 
pressures larger than zero. The two-dimensional Gaus- 
sian core model, which approaches the hard disk system 
in the limit of low temperature and low density [30j , has 
been used previously to study the melting transition in 
two dimensions [30] ISlj. 

To make contact between numerical calculations in 
the particle system and continuum elasticity theory, we 



determine, at T = 0, the displacement field [32] 

u(r,) = v[ - r, (2) 

caused by the introduction of the defect into the system. 
Here, y[ and denote the position of particle i with 
and without the defect, respectively. The displacemet 
field completely describes the response of the system's 
structure to the perturbation introduced by the defect. 
Numerically, we determine displacement fields by insert- 
ing a particle into or removing it from a perfect crystal 
on a triangular lattice. The system is then relaxed to a 
new minimum energy configuration by steepest descent 
minimization at constant volume of the simulation box. 
Periodic boundary conditions apply. Typically, about 
tens of thousands of steepest descent steps are required 
to determine minimum energy structures accurately. In 
each minimization step, each particle is moved in the 
direction of the force acting on the particle where the 
absolute value of the displacement in chosen to be small 
enough to ensure that the energy of the system is re- 
duced in each step. The displacement u(ri) of particle i 
is then simply the vector which connects the position of 
particle i before the minimization, r^, with its position 
after the minimization, v[. 

The largest system we study here consists of = 
199, 680 Gaussian core particles (without the extra par- 
ticle) at a number density oi p — 0.6(t~^ correspond- 
ing to a lattice constant of a = 1.3872cr. The al- 
most square simulation box has length Lx = 416a and 
height Ly = (\/3/2)480a = 415.692a with aspect ratio 
Ly/Lx = 0.99926. 

2.1 Constant V or constant pi 

In calculating the displacement fields caused by point 
defects the question naturally arise whether one should 
do that at constant volume V or at constant pressure p. 
Naturally, the choice should depend on the particular 
experimental situation one is interested in. As we will 
show here, however, the displacement fields caused by a 
point defect at constant pressure and at constant volume 
are simply related. To determine how they are related, 
consider a perfect triangular crystal at T = enclosed in 
a rectangular cell of volume Vq with appropriate aspect 
ratio. For this particular volume, the crystal is under the 
hydrostatic pressure po- Insertion of a point defect into 
the crystal at a fixed total volume distorts the crystal 
and atom i is displaced by 

uy„(r,)=r:(yo)-r,(yo), (3) 

where the subscript Vq indicates that the displacement 
field \iVa{Yi) is obtained at constant volume Vb- In the 
above equation. r';(Vb) and ri(Vo) are the positions of 
atom i in the system of volume Vq with and without 
the defect, respectively. If one requires, however, that 
the defects is created at constant pressure poj the vol- 
ume of the simulation cell changes from Vq to Vi (typi- 
cally, it will increase for an interstitial and decrease for 
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a vacancy) and the atoms are displaced by a different 
amount, 

Up,(rO=r',(Fi)-r,(Fo), (4) 

where the subscript po implies that the displacement 
field is considered at constant pressure. Note that here 
we assume that during the generation of the defect the 
simulations cell only expands or contracts, but does not 
change its shape. This assumption can be lifted as dis- 
cussed below. We now imagine that the defect gener- 
ation at constant pressure is carried out in two steps: 
first the system is homogeneously dilated without de- 
fect from volume Vq to volume Vi; in the second step, 
the defect is inserted at constant volume Vi. This two 
step operation corresponds to adding and subtracting 
ri(Vi), i.e., the position of atom i at volume Vi in the 
absence of the defect, to the right hand side of the above 
equation, 

Up„(r,) = r,(Fi) - r,(^i) + r',{V,) - r,(Fo). (5) 

What one obtains in this way is 

Upo ) = uy, (r, ) u,, (r, , Vo , ) , (6) 

where uvi(ri) = r-(Vi) — ri(Vi) is the displacement field 
obtained by inserting the defects at volume Vi for fixed 
simulation cell and u/i (r^ , Vq , Vi ) = ri(Vi) — ri(Vb) is the 
displacement field corresponding to a homogeneous di- 
latation (or contraction) of the perfect crystal without 
defect from volume Vq to volume Vi . This simple defor- 
mation corresponds to a displacement Uh{ri,Vo,Vi) = 
{Vi/Voy^^Ti (in three dimensions the exponent is 1/3). 
Hence the displacement fields for constant pressure and 
constant volume are related by: 

[Vi 

Upo(r,) = uv,(ri) -I- W— r;. (7) 
V Vo 

Thus, one can determine the constant-pressure displace- 
ment at pressure po by calculating the constant-volume 
displacement at volume Vi , the volume at pressure po in 
the presence of the defect. 

Similar considerations can be used to relate the 
constant-pressure and constant-volume displacement 
fields if the simulation cell is permitted to change shape 
as well as volume during the constant-pressure defect in- 
sertion. In this case, the simulation cell is characterized 
by the vectors a and b along its edges [33] . It then turns 
out that the displacement field of a defect inserted into 
an initially rectangular simulation cell with edge vec- 
tors a and b at constant hydrostatic pressure is simply 
related to the displacement field for fixed cell vectors 
a' and b' which, in general, differ from a and b. For 
sufficiently large systems, however, a fixed shape of the 
simulation cell is only a very weak constraint. In par- 
ticular, a displacement field which tends to be isotropic 
at large distances may lead to a change in aspect ratio 
of the simulation cell at constant pressure, but not to 
a change in the relative orientation of the edge vectors. 
All calculations of this paper arc carried out for fixed 
and nearly square simulation cells. 



2.2 Interstitials 

We first determine the displacement field of a single in- 
terstitial particle. This type of point defect can exist 
in different configurations [4] with displacement fields 
of different symmetries [53]. The three lowest energy 
structures are shown in Fig. [TJ In one minimum-energy 
configuration, termed I2 interstitial and shown in Fig. 
[1^, the extra particle and one of the original particles 
arrange themselves at equal distance around the lattice 
position of the original particle leading to a two-fold 
symmetry. This is the two dimensional analogue of the 
crowdion in an fee crystal [35]. The displacements are 
largest on the main defect axis, which can be aligned in 
any of the three low-index directions of the lattice. Since 
the defect symmetry differs from that of the underlying 
triangular lattice, one may wonder whether the rectan- 
gular periodic boundary conditions used in the calcula- 
tion favor the two-fold defect symmetry. Calculations 
carried out with hexagonal boundary conditions, how- 
ever, yield identical results demonstrating that the de- 
fect symmetry is not imposed by the symmetry of the 
boundary conditions. 

Other low-energy defect configurations include the 
I3 interstitial with three-fold symmetry shown in Fig. 
[Dd and the Id interstitial or dumbbell interstitial shown 
in Fig. (TJ;. In the I3 configuration, the extra particle is 
located at the center of a triangle spanned by three near- 
est neighbor lattice points and the surrounding particles 
are displaced outward with respect to their original po- 
sitions. In the dumbbell configuration, the interstitial 
particle and one of the original particles compete for 
one lattice position as in the I2 interstitial, but the line 
connecting them is orthogonal to one of the low-index 
lattice directions. In contrast to the I2 interstitial, the 
Id is not concentrated on one single axis. 

The interstitial configurations observed in the Gaus- 
sian core model have energies that differ by less than 
0.1% of the total defect energy. These energy difference 
correspond to roughly 20% of the thermal energy ksT at 
melting. At finite temperatures that are not too low, in- 
terconversion between the various defect configurations 
is facile and all three of them play an important role 
during defect diffusion [16] . 

2.3 Vacancies 

Also vacancies can occur in various configurations with 
displacement fields displaying quite complex patterns 
and symmetries lower than that of the underlying lat- 
tice. Three minimum energy configurations are shown 
in Fig. [2] In the the vacancy configuration V2 (see Fig. 
[2^), particles move mainly on the a;-axis to partially fill 
the void left by a removed particle. As a result, particles 
above and below the void site move outward generating 
a vortex-like displacement field. This two-fold vacancy 
V2 has the same symmetry as the I2 interstitial, but 
their displacement patterns are not simply related by 
inversion. In particular, the vortex structure observed 
for the vacancy is absent in the interstitial case. In the 
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Figure 1: Displacement fields (left hand side) and par- 
ticle configurations (right hand side) of the I2 intersti- 
tial (a), the /a interstitial (b) and the Id interstitial (c). 
The arrows representing the particle displacements are 
exaggerated in length by a factor of 20 for better visibil- 
ity. On the right hand side, the small grey spheres indi- 
cate the sites of the perfect triangular lattice. The blue 
spheres represent particles with 6 neighbors according to 
the Voronoi construction (black lines). Yellow spheres 
are particles with 4 neighbors, orange and green spheres 
represent particles with 5 and 7 neighbors, respectively. 



configuration V3 with threefold symmetry, particles par- 
tially fill the vacancy void by moving in along three axes 
rather than two. On the other three low-index axes, 
particles are moved outward in response to the removed 
particle. A vacancy configuration analogous to the Id 
interstitial seems not to be stable even at T = 0. A 
configuration prepared in this symmetry ends in an an- 
tisymmetric configuration VJj (see Fig 121;). Energetically, 
configurations V2 and Va are equal and lower than config- 
uration V3 by more than twice the thermal energy k^T 
at melting thus exceeding the energy difference of the 
corresponding interstitial configurations by more than 
an order of magnitude. This energy difference is less 
than 10% of the total defect energy. 



Figure 2: Displacement fields (left hand side) and par- 
ticle configurations (right hand side) of the V2 vacancy 
(a), and the V3 vacancy (b) and the Va vacancy (c). The 
arrows representing the particle displacements are exag- 
gerated in length by a factor of 20 for better visibility. 
The color code is the same as in Fig[TJ 

3 Elasticity Theory 

Near the defect site non-linearities and discrete lattice 
effects dominate the displacement pattern as evidenced 
by the highly anisotropic local structure of vacancies and 
interstitials. Far away from the defect, however, the 
perturbation of the 2d-crystal should be described accu- 
rately by continuum elasticity theory. In this regime, the 
response of the system to a point defect should depend 
on the specific form of the interaction potential only 
through the particular values of the elastic constants. 
To verify to which extent elasticity theory is valid for 
two-dimensional colloidal crystals of soft particles, we 
first review the basic equations of elasticity theory and 
then solve them for an idealized singular defect model 
consisting of a pair of singular forces of equal magnitude 
and opposite direction [36l [371 [38] . 

Linear elasticity theory is usually formulated in terms 
of the symmetric strain tensor [52] 
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where Ui and are the i-ih component of the displace- 
ment and the position, respectively. For small strains, 
Hook's law applies and the stress cr,;j is linearly related 
to the strain e- 



ijkl ^kl ; 



(9) 



where Cijki is the stiffness tensor. Here and in the 
following, summation over repeated indices is implied. 
For isotropic materials, such as two-dimensional crys- 
tals with triangular lattice, this relation reduces to 



(10) 



where A and fi are the so-called Lame coefficients. The 
Lame coefficient is also called the shear modulus. 

In order to calculate the displacement field generated 
by a point defect one must be able to determine how the 
elastic continuum reacts to external forces. The condi- 
tion that the forces on each infinitesimal volume element 
balance leads to 



9r, 



+ = 0, 



(11) 



where fi is component z of a given volume force /(r) 
acting at r. Using the stress-strain relation from Equ. 
(|f op . these equilibrium conditions can be formulated in 
terms of the strains rather than the stresses, 



X^ekk + + /, = 0, 



(12) 



Inserting the definition of the strain into this equation 
one obtains the equilibrium conditions for the displace- 
ment field u(r). 



(A 4- /i) ^ ^"-^ 



nAui + fi = 0, 



(13) 



where A = d^/dx^ + d^/dy^ is the Laplace operator. 
Solving this equation for a particular arrangement of 
forces used to model the point defect then yields the 
displacement field caused by the forces. 

From the displacement field one can then determine 
the energetics of the point defect. In terms of the strain 
tensor and the Lame coefficients the elastic free energy 
density of the system is given by 



A 



9 



2 I 2 



2 ' ij 

Accordingly, the energy density at T = is given by 
A 



(14) 



(15) 



The last term of this equation stems from the work done 
against the pressure p by the dilatation e^k ■ The strain 
tensor can also be written as the sum of a trace-free shear 
and a homogeneous dilation leading to the expression 



SijCkk 



K 



(16) 



where K is the so-called bulk modulus related to A and 
M by 

A' = A + fi. (17) 

The Poisson ratio f, i.e., the negative ratio of transverse 
strain to axial strain upon uniaxial loading, is given by 



A 



X + 2fi^ K + n 



(18) 



and describes how a material reacts when stretched. In 
the next subsection we will calculate the elastic con- 
stants for our system at T = 0. 

3.1 Elastic moduli 

For a crystal in which particles interact with a pair po- 
tential v{r) depending only on the interparticle distance 
r the total energy E oi N particles is given by 



E 



r,|), 



(19) 



where and are the positions of particles i and j 
respectively. In this case and for T = 0, the energy 
density eo of the undistorted lattice, the pressure p, as 
well as the clastic constants K and ^ can be calculated 
from simple lattice sums: 



eo 



-y 

2 ^ 



and 



IE' 



\ 2 4 

^ (^0 ( — ) (^0^ 



(20) 



(21) 



(22) 



(23) 



Here, v'{r) and v"{r) are the first and second deriva- 
tive of the pair potential, respectively, p is the number 
density, is the distance of particle i from the origin, 
and Xi and yi are its Cartesian coordinates. The lattice 
position is chosen such that there is one particle at the 
origin. The sums in the above equations must include 
a sufficient number of particles to ensure convergence of 
the sums. (The prime on the sum symbol indicates that 
the particle at the origin is not included in the sum). 

The elastic constants ^, and K as well as the Pois- 
son ratio the pressure p and the energy density e, 
calculated using such sums, are shown in Fig. [3] as a 
function of the lattice constant a, which, in a triangu- 
lar lattice, is related to the density by p = 2/^/30^. At 
a density of p = O.Qa^, the density at which all cal- 
culations presented in this article are carried out, the 
elastic constants have the values K = 1.2089 ecr"^ and 
^ = 0.060183 £0-^2, the pressure is p = 0.544245 ecr^^, 
and the Poisson ratio is = 0.905151. The energy den- 
sity is eo = 0.269125 ecr"^ corresponding to an energy 



5 



1.2 
1 

0.8 
0.6 
0.4 
0.2 




1\ \ 
— \ 1 


0.08 




/\ 

\ 

\ 

\ 


1 \ 

~ \ \ 


\ 0.06 
\ 0.04 
A 0.02 

\ -. 

\ \ 


- / 
" / 

- l' 


— \ \ 


1 


2 3 
a 




\ \ 








\ \ \ 






1 1 1 


1 T-~-tn 




ad_J 1 1 1 1 



2 3 
a 



Figure 3: Bulk modulus K (solid line), shear modulus 
fi (dashed line), pressure p (dash-dotted line), energy 
density e (thin solid line), and Poisson ratio i' (dotted 
line) as a function of the lattice constant a. The moduli 
are given in units of e/cr^ and the lattice constant in 
units of a. In the inset, the shear modulus fi is displayed 
on a larger scale. The vertical thin dotted line indicates 
the lattice constant a ~ 1.3872cr corresponding to the 
density p = O.Gcr^^, at which all calculations discussed 
in this article are carried out. 

per particle oi E/N = 0.448542e. Note that at this 
density, the system is stabilized against shear only by 
interactions beyond nearest-neighbor contributions; es- 
timation of /X from nearest neighbor interactions only 
yields a negative shear rate at this density. While the 
bulk modulus K increases monotonically with the den- 
sity (and, as the pressure p, is proportional to for 
small densities), the shear modulus reaches a maximum 
at a 1.67 and then rapidly decays to zero for lattice 
constants larger and smaller than that. This behavior 
of the shear modulus is a reflection of the phenomenon 
of reentrant melting observed in the three-dimensional 
Gaussian core model [26t [27] and indicates that also in 
two dimensions a Gaussian core crystal melts if suffi- 
ciently compressed. 

3.2 Point defect model 

We next use elasticity theory to determine the displace- 
ment field created by introducing an idealized point de- 
fect into a perfect isotropic crystal. The dilatation (or 
contraction) caused by the defect is modeled by two or- 
thogonal pairs of forces. Each pair consists of two forces 
of equal magnitude F but opposite directions acting at 
two points separated by the distance h. If one assumes 
that one force pair acts in x-direction and the other one 
in y-direction and that the defect is centered a the origin, 
the total force density is given by 



f(r) = -F6{r)e., 
-FS{r)ey 



FS{r - he^)e.j, 
FS{r — hey)ey. 



(24) 



Here, 6; 

direction, respectively, and 6{r) is the Dirac delta- 



and By are the unit vectors in x- and y- 



function in two dimensions. One then lets the separation 
h go to zero and the force F go to infinity in a way such 
that Fh remains constant. This defect model, in which 
the net force acting on the material vanishes, is equiva- 
lent to inserting a small circular inclusion into a hole of 
different size [37] . 

The displacement field caused by this type of point 
defect can be determined by first calculating the Green's 
function for a singular force and than carrying out the 
limit h ^ 0. Alternatively, one can carry out the limit 
h ^ first and then solve the equilibrium condition for 
the force density obtained in that way. In the follow- 
ing we will calculate the displacement field of the point 
defect model using this second approach, in which pe- 
riodic boundary conditions can be taken into account 
particularly easily. 

Carrying out the limit /i ^ as described above the 
total force density of Equ. ((24]) reduces to 



f(r) = -FhVS(r). 



(25) 



Inserting this expression into Equ. (|13p one obtains the 
equilibrium condition for this simple point defect model. 



ori or^ 



pAui 



Taking the divergence on both sides yields 
A{\ + 2p)-^ = FhA6{r). 



(26) 



(27) 



To solve this equation it suffices to find a displacement 
field that obeys 



{X + 2p)^^Fh5{r) 



(28) 



Using the Helmholtz-decomposition in two dimensions, 
we now write the displacement in terms of the gradients 
of two scalar functions (^(r) and A{r) as a sum of an 
irrotational and a divergence-free part, 



dA 



dri 



(29) 



where the matrix ojij exchanges the components of the 
gradient and changes the sign of one of them: ujn = 
UJ22 = and —UJ21 ~ ^12 = 1- Then, Equ. (|28p becomes 



A0(r) = 2Trj6{r), (30) 

where we have used the fact that ojijdA/drj is 
divergence- free and the parameter 7, which has the di- 
mension of an area, is given by 



Fh 



1 



2ti{\ + 2p) 



(31) 



Equation (|5(I]) is the Poisson equation of electrostatics 
for a point charge of strength —7 in two dimensions. 

A similar equation can be derived for the the scalar 
function A{r) by taking the 2d-vorticity, defined as 



6 



ujijdvj / dri for an arbitrary vector field v — {vl^V2)^ of 
both sides of Equ. (|26p . Since the vorticity of a gradient 
field vanishes, one obtains the biharmonic equation 



^A(A^(r)) = 0. 



(32) 



This equation is obeyed if the scalar field A{v) is a so- 
lution of the Laplace equation 



AA(r) 



0. 



(33) 



In the following, we will use the trivial solution A{v) = 
const and satisfy the boundary conditions through 
proper solution of the Poisson equation ([50]) for the 
scalar field (j>. 

To do that, we note that K{i-) ~ ln(r)/27r is a so- 
lution of AK = 6{r) (see, for instance, Ref. [3S])j and 
hence we obtain the Green's function 



(r) = 7 ln(r). 



(34) 



The corresponding displacement field u(r) follows by dif- 
ferentiation according to Equ. ([M]) , 



(35) 



Thus, the displacement field caused by the point de- 
fect is isotropic and long-range with a magnitude that is 
proportional to 1/r. This result is valid for an infinitely 
extended elastic medium where the boundary conditions 
u = apply at infinity. This situation, however, does 
not correspond to the boundary conditions applied in 
computer simulations. In the following section we will 
discuss how to solve Equ. (j30p with the appropriate 
boundary conditions. 

3.3 Boundary conditions 

In comparing the results of particle simulations with 
those of elasticity theory it is important to realize that 
the displacement fields predicted by continuum theory 
are of long-range nature. Therefore, it is crucial that cor- 
responding boundary conditions are used in both cases. 
All simulations discussed in this paper are done with 
periodic boundary conditions in order to minimize finite 
size effects and preserve the translational invariance of 
the perfect lattice. Hence, also the continuum calcu- 
lations need to be carried out with periodic boundary 
conditions. For a rectangular system with side lengths 
Lx and L^, periodic boundary conditions require that 
u(r) = u(r-l-l), where 1 = {iLx,jLy) is an arbitrary lat- 
tice vector with integer i and j. In the following, we will 
solve the Poisson equation ([50]) for this type of boundary 
conditions. 

We start by noting that the homogenous part of the 
Poisson equation (|30p admits the non-trivial solution 
(/)o(r) = const that satisfies the boundary conditions. 
Therefore, one needs to consider the extended Green's 
function for the solution of the general Poisson equation 
A0(r) = 27r/9(r) [3S1I1D]. In this case, the right hand 



side of the Poisson equation must be orthogonal to the 
solution (j)o{r), 



dr (^o(r)p(r) = const / drp(r 



0. 



(36) 



In electrostatics, this condition corresponds to charge 
neutrality (the physical meaning of this condition in our 
case will be discussed below). To satisfy this orthogo- 
nality condition we must modify the Poisson equation 
([50]) by subtracting 1/A from the delta function. 



A0(r) = 27r7 



Sir) 



(37) 



where A is the area of the rectangular basic cell. In this 
modified equation, the right hand side contains a ho- 
mogeneous "neutralizing background" that exactly com- 
pensates for the "charge" of the delta function. Solution 
of this equation yields the extended Green's function of 
the problem. To obtain a unique solution (/)(r) of this 
equation one must furthermore require that this solu- 
tion be orthogonal to 0o(r)j 



dr <po{r)(f){r) ^ const / dr0(r) = 0. (38) 



For our case this condition is irrelevant, as only deriva- 
tives of (j){r) carry physical significance. Once the func- 
tion (f>{i-) has been determined by solving Equ. ([57]) . the 
displacement field follows by differentiation. 

3.3.1 Rigid circular container 

Before we embark on the solution of the extended Pois- 
son equation (|37p for periodic boundary conditions, we 
illustrate the concepts introduced above by determin- 
ing the displacement field of a point defect in an elastic 
material enclosed in a container with rigid walls. Due 
to these walls, the component of the displacement field 
normal to walls must vanish at the wall, u±_ = 0. No 
condition applies for the parallel component For a 
rectangular container such rigid wall boundary condi- 
tions are equivalent to periodic boundary conditions. If 
we assume, without loss of generality, that the point de- 
fect is located at the center of the rectangular periodic 
cell, the component of the displacement field normal to 
the boundary of the periodic cell vanishes also in this 
case. In the following, we will determine the effect of 
such rigid boundary conditions on the displacement field 
of a point defect located at the center of a circular cav- 
ity enclosed by hard walls. For this case, which exhibits 
all complications mentioned above, a simple analytical 
solution can be easily obtained. 

Consider a two-dimensional elastic isotropic material 
enclosed in a circular container of radius R. We choose 
the coordinate system such that the origin is at the cen- 
ter of the container. To determine the displacement field 
caused by a point defect of strength 7 placed at the 
origin we need to solve Equ. (j37p under the condition 
that the at distance R from the origin the displacement 
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u± normal to the wall vanishes. We construct a solu- continuum calculations agree very well for all defect dis- 

tion by superposing the solution for the extended ma- tances larger than about 15 lattice constants. In Fig. 

terial, Uoo(r) ~ 7r/r^, and a homogeneous contraction. [5] we depict the relative error which we define in the 

Uc(r) = — ar, following way: 



u(r) 



(39) 



While Uc(r) corresponds to a homogeneous contraction 
without shear, Uoo(r) corresponds to a pure shear with- 
out dilatation (except for r = 0). To satisfy the bound- 
ary conditions at r = i?, we set a = j/R^ obtaining 



u(r) 



(40) 



This displacement field corresponds to the "potential" 



(r) = 0oo(r) + Mr) = ^ i^^^Jl - Ij^ + I 



(41) 



where the last constant on the right hand side takes care 
of the condition expressed in Equ. p8|) . It is straight- 
forward to verify that 

A<^(r) = AM (r) + A0c (r) = 2TT^S{r) - ^ , (42) 



such that the potential from Equ. (|4T|) satisfies the ex- 
tended Poisson equation ([57]) . 



Simulation 

Continuum Theory 

y/r 




^ |up(r) -Uc(r)| 
" |uc(r)| ■ 



(43) 



Here, Up(r) and Uc(r) are the displacement fields ob- 
tained from the particle system and from the predictions 
of continuum theory, respectively. The relative error 
close to the defect can be larger than 100% (red). For 
distances of r > 20a with find an error of approximately 
1 — 5% (green). Close to the rigid container the relative 
error increases again due to discrete lattice effects. 
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Figure 5: Color coded relative deviation ^ calculated 
according to Equ. (|43p of an I2 interstitial at the origin 
of a rigid circular box with radius R = 108.19cr. The fit 
parameter 7 = 0.234495(7^ was found by minimizing the 
sum of the relative error of particles at distances larger 
than 30(7. The contour lines in white, black, and blue 
represent an error of 1%, 2%, and 5%, respectively. 



Figure 4: Displacement components Ux and func- 
tion of the distance r (thin solid lines) for an I2 inter- 
stitial with its main axis oriented in a;-direction. The 
circular hard wall container has a radius of i? = 106.8(7. 
Also plotted is the displacement computed from contin- 
uum theory according to Equ. ([3D]) (dashed line) for a 
defect strength of 7 = 0.234495cr^ which best fits the 
numerical results in the far field and the simple 1/r- 
behavior (dotted line). 

A comparison of results obtained numerically for an 
I2 interstitial and the prediction of elasticity theory 
(Equ. (|40|) ) is shown in Fig. |4l For the particle sys- 
tem the rigid container was realized by carrying out the 
calculation in a larger system in which all particles be- 
yond a distance of R from the origin where kept at fixed 
positions. The displacements obtained from particle and 



Similar agreement is found also for the energy density 
as shown in Fig. \6\ For the displacement field of Equ. 
(|40|) one finds, using Equ. p5p the energy density 



e(r) 



7' 



2K^ 

i?4 



(44) 



This prediction of elasticity theory matches the energy 
density determined numerically in y-direction (see Fig. 
|6|). Due to the strong anisotropy of the I2 defect, larger 
deviations are observed in cc-direction. For distances of 
more than about 15 lattice spacings the energy density 
plateaus at a constant value. In this regime, the en- 
ergy density is essentially constant, e = 2pj/Rl^, and 
corresponds to the work done by the defect against the 
pressure p. As discussed below, the plateau value of the 
density is related to the neutralizing background on the 
right hand side of Equ. ([57]) . 
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Figure 6: Energy density e as a function of distance 
r measured along the x-axis (solid line) and the j/-axis 
(dotted line) for the particle system in a circular con- 
tainer with hard walls and radius R ~ 71.91(7. Also 
shown is the energy density calculated from continuum 
theory according to Equ. (dashed line) . 



3.3.2 Ewald summation 

For an isolated defect at the center of a rectangular cell, 
the symmetry imposed by periodic boundary conditions 
requires that the components of the displacement field 
orthogonal to the surface of the cell vanish. In a first at- 
tempt to obtain the displacement field for such boundary 
conditions one may start from the solution for the ex- 
tended material and satisfy the boundary conditions by 
placing "image defects", each of which carries the dis- 
placement field for the infinitely extended material, at 
appropriate positions. For a rectangular cell, an infinite 
number of image defects arranged on a regular lattice 
with lattice constants and Ly in x- and y-direction, 
respectively, are required. These image defects, which 
are analogous to the image charges of electrostatics, cor- 
respond to the defects in the periodic images of the basic 
simulation cell. Superposition of the displacement fields 
of all image defects then yields the displacement field for 
periodic boundary conditions. 

Due to the long-range nature of the defect field for 
the infinite material, however, such a summation of the 
contribution of all image defects leads to displacement 
fields that are only conditionally convergent. A more 
appropriate treatment that avoids this problem consists 
in determining the Green's function of the Poisson equa- 
tion (|37p for periodic boundary conditions. The require- 
ment imposed by the periodic boundary conditions can 
be easily satisfied by expressing the solution as a Fourier 
series and solving the Poisson equation in Fourier space. 
This treatment, however, leads to series that are only 
conditionally convergent with values that depend on 
the summation order. The solution of this problem us- 
ing so called Ewald sums is known from electrostatics 
[4T1 l42l [40] and consists in separating the conditionally 
convergent series into a real space and and a Fourier 



2tt X ^ 



fc2 



Ir + l|- 



■ cos(k • r) 



(45) 



Here, Ei{x) = J^^{e* /t) dt is the exponential integral 
and A is the area of the rectangular cell. The first sum is 
over all lattice vectors 1 and the second sum is over all re- 
ciprocal vectors k consistent with the periodic boundary 
conditions. The adjustable parameter rj^ set to a value 
of T] = 6/Lx here, determines the rate of convergence 
of the two sums, but the value of the sums does not de- 
pend on rj. The exclusion of the k = term in the above 
equation stems from the requirement that both 0(r) and 
the right hand side of the Poisson equation need to be 
orthogonal to (/'o(r) as expressed in Eqs. ((36)) and p8|) . 
It is easy to show by direct calculation of the Lapla- 
cian At/) that the above expression for (j){r) indeed obeys 
Equ. ([57|) and thus implies a neutralizing background 
of magnitude 7/A as in the previous example. 

From Equ. ([45]) for the scalar function (f){r) the dis- 
placement field of a point defect in a system with peri- 
odic boundary conditions is determined by differentia- 
tion, 



|r+l|- 



27r 



A ^ fc2 
k#0 



-ki sin(k • r) 



(46) 



For the systems and parameters considered in this paper, 
the real space sum can be truncated after the first term 
and the Fourier space sum can be evaluated accurately 
using about 50x50 reciprocal vectors. 

As mentioned above, exactly the same boundary con- 
ditions apply if the system is constrained by hard walls 
to reside in an area of given size. Also in that case, 
the boundary conditions require that the component of 
the displacement field orthogonal to the walls vanishes. 
Hence, in the continuum description, hard walls have the 
same effect as an infinite array of image "charges" (plus 
"neutralizing background") placed on a regular lattice 
with a geometry determined by the wall positions. The 
effect of such image charges and the neutralizing back- 
ground is, therefore, not a pure artifact of the periodic 
boundary condition applied in the simulations, but oc- 
curs also in experimental realizations of colloidal crystals 
of purely repulsive particles which need to be kept to- 
gether by confining walls. Accordingly, the analysis of 
displacement patterns (and defect interactions) observed 
experimentally requires a similar treatment as that used 
here for the interpretation of our simulation results. 
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3.4 Electrostatic analogy 

In electrostatics, the technique of Ewald summation is 
used to determine the energetics of periodic systems con- 
taining point charges and dipoles. When one uses this 
technique, one imphcitly stipulates that the charges are 
immersed in a homogeneous background that compen- 
sates for the point charges and establishes overall charge 
neutrality. This neutralizing background is imposed by 
the periodic boundary conditions; without it, no peri- 
odic solution of the Poisson equation exists. Since math- 
ematically the situation we face when determining the 
displacement field of point defects is identical to that of 
electrostatics, one may wonder about the physical mean- 
ing of the neutralizing background in our Equ. (|37|) . 

To address this question, we note that the local vol- 
ume change, or dilatation, due to a displacement field 
u(r) is given be the trace e^^ of the corresponding strain 
tensor [32]. The total change in volume AV of a certain 
region G is then given as the integral over the dilatation, 

AV= [ drefefe(r). (47) 
Jg 

On the other hand, it follows from the definition of (j> 
(see Equ. (HH])) that the trace of the strain tensor is 
equal to the Laplacian of 0, 

A0(r) - efcfc(r). (48) 

Thus, the Poisson equation ([57]) is an equation for the 
local dilatation. According to this equation, at the de- 
fect site the dilatation is required to have a delta like 
peak of strength 2tt^. The total volume change caused 
by this singular dilatation, AV = dr 27r7(5(r) = 27r7, 
is exactly compensated by the total volume change 
originating from the constant neutralizing background, 
AV = — J^dr2Tr-f/A = —I-k^. (Calculating the strain 
tensor directly from the displacement field of Equ. (j46ll 
indeed yields e^fc ~ —2ti^ /A.) Therefore, the condition 
of charge neutrality of electrostatics corresponds to the 
requirement of constant volume in our case. In this anal- 
ogy, the charge density of electrostatics corresponds to 
the local dilatation (the charge corresponds to the vol- 
ume change) and the role of the electric field is played 
here by the displacement field. 

This interpretation of the Poisson equation ([571) ^l^o 
suggests a definition of the defect volume Vd- As men- 
tioned above, introduction of a point defect of strength 
7 leads to a total volume change which can be viewed 
as the volume of the defect, 

Vd = 2^7. (49) 

With periodic (or rigid) boundary conditions the system 
as a whole is prevented from changing volume and the 
volume change due to the defect is exactly compensated 
by the homogeneous neutralizing background. This de- 
fect volume is also what one gets when calculating the 
expansion of a circle under the 1/r-deformation caused 
by the idealized defect model (see Equ. (|35p). Measur- 
ing the parameter 7, for instance by fitting the displace- 
ment field far from the defect to the continuum theory 



results, thus permits to determine the defect volume. 
For the I2 at a density oi p = 0.6cr~^, for example, we 
found a volume of Vd = 1.44(T^, which is slightly smaller 
then Vq = 1.66(7^, the volume per particle in the perfect 
lattice. 

The neutralizing background appearing in the Pois- 
son equation (jST]) also figures in the energy density and 
can contribute considerably to the total defect energy. 
According to Equ. ([15]), the energy density includes the 
term Cp = —^kkP arising from the work carried out by 
the defect against the pressure p. Away from the singu- 
larity at the origin, this component of the energy density 
is constant, Cp = 2TT'y/A, and it dominates for large dis- 
tances from the defect. Although Cp is proportional to 
1/A and therefore small in general, integrating it over 
the entire area (leaving out the unphysical singularity 
at the origin) yields an energy contribution of Ep = 2717 
which is independent of system size and can be substan- 
tial. For an I2 interstitial in our Gaussian core system at 
p ~ 0.6a~^, for instance, this contribution amounts to 
more than 50% of the total defect energy. Interestingly, 
no such pressure contribution arises for the displacement 
field of Equ. ([35]) obtained for the infinitely extended 
material. Thus, the condition of fixed volume imposed 
by the periodic boundary conditions leads to measurable 
effects also in the large system limit and even boundary 
conditions applied at infinity matter. 

4 Comparison of simulation and 
continuum theory 

In this section we compare the results of the particle- 
based simulations with the predictions of continuum 
elasticity theory obtained in the previous section. In 
particular, we verify at which distances from the point 
defect elasticity theory becomes valid and which effect 
boundary conditions have on the displacement fields. 
We first consider the displacement fields of interstitials, 
then those of vacancies. 

4.1 Interstitials 

As discussed in Sec. [2l insertion of an interstitial particle 
into a perfect lattice can lead to different displacement 
patterns, all of which are highly anisotropic near the de- 
fect site. Farther away from the defect the anisotropy 
should subside as the isotropic behavior expected from 
elasticity theory sets in. This is indeed what is observed 
for intermediate distances from the defect as shown in 
Fig. [7] for an I2 defect. In the bottom panel of this 
figure, the displacement magnitude |u(r)| is plotted as 
a function of the distance r from the defect. Each dot 
corresponds to one particular particle. While for short 
distances the displacement magnitude is not a unique 
function of r due to the anisotropy of the defect, at larger 
distances |u(r)| is essentially determined by r. In this 
intermediate regime, the displacement magnitude seems 
to follow the 1 /r-behavior predicted by elasticity theory 
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for the infinitely extended material. Due to the periodic 
boundary conditions, however, the displacement mag- 
nitude cannot remain isotropic as the boundary is ap- 
proached. In fact, the periodic boundary conditions lead 
to a spread of |u(r)| at larger distances. The two prongs 
observed in the bottom panel of Fig. [7] correspond to 
the directions along the x- and y-axes and along the 
diagonals. 




r 



Figure 7: Top: Angle 9 between the displacement vector 
u and the position vector r as a function of the distance 
r from the defect for an I2 interstitial in a 208x240 
particle system. Each dot corresponds to one particle. 
Bottom: Displacement magnitude |u(r)| as a function 
of r. Also shown as a dashed line is the 7/r-line for 
7 = 0.23165cr2. 

This kind of behavior is observed even more clearly 
for the displacement directions. In the top panel of Fig. 
[71 the angle 9 between the displacement u(r) and the po- 
sition vector r is plotted as a function of distance r. As 
in the bottom panel, each dot corresponds to an individ- 
ual particle. For an isotropic displacement field, the dis- 
placement and the position vector are perfectly aligned 
and 9 = 0. Thus, non-zero angles 6 are an indication of 
anisotropy. For small distances r angles 6 larger than 
7r/4 occur. For intermediate distances, 20 < r < 150, 
the angle 9 is small since in this regime the displace- 
ments approximately points away from the origin. For 
larger distances, the periodic boundary conditions then 
lead to a spread in 9 and deviations of up to = 7r/2 
are possible. 

The long-distance behavior described above is per- 
fectly reproduced by linear elasticity theory. As shown 
in Fig. [5] for I2 interstitials and in Fig. [3] for I3 and 
Id interstitials, respectively, the displacement calculated 
using Ewald summation according to Equ. (|46p agrees 
very well with the numerical results for all distances 
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Figure 8: Displacement components and Uy of the 
I2 interstitial as a function of distance r along the x- 
axis and y-axis, respectively (solid lines). The main axis 
of the I2 defect is oriented in x-direction. Also plotted 
are the displacement computed from continuum theory 
by Ewald summation according to Equ. (|46l) (dashed 
line), and simple 1/r-behavior (dotted line). For short 
distances, the behavior of the displacement along the 
X-axis is exponential (dash-dotted line) as described by 
a simple bead-spring model [34j . The inset shows the 
region close to the defect location. A defect strength of 
7 = 0.23165(7^ was used here since this value yields the 
best fit of the results obtained from elasticity theory and 
the numerical results at large distances from the defect. 

larger than about 10 — 15 lattice spacings. In par- 
ticular, the deviations from the 1/r-behavior near the 
cell boundary arc perfectly captured by elasticity the- 
ory with periodic boundary conditions. For small dis- 
tances, on the other hand, the displacement field in the 
particle system is highly anisotropic with strong devia- 
tions between the x- and y-direction. In this non-linear 
core region elasticity theory is not applicable and the 
displacements of the three configurations differ. The ex- 
ponential short-range dependence of on the distance 
r for I2 interstitials is, however, captured by a simple 
bead-spring model discussed in Ref. |34| . In this model, 
the exponential decay constant can be related to the 
elastic constants of the material. 

In the comparison of the results obtained for the par- 
ticle system with those of continuum theory the defect 
strength 7 is treated as an adjustable parameter. For 
each configuration, I2, I3, and Id, the particular dis- 
placement strength 7 was found by optimizing the rela- 
tive error (see Equ. at distances larger than 30. Oct 
from the origin of the defect. For the I2 defect, a value of 
7 = 0. 23165(7^ yields the best fit. I3 and Id interstitials 
produce a slightly larger displacement with a strength of 
7 = 0.23470-2 and 7 = 0.24250-^, respectively. The ques- 
tion arises if this fit is independent on the size of the 
box. In Fig. [10] we depict the displacement of an I2 de- 
fect for different box sizes together with the results from 
the Ewald summation. The displacements plotted in the 
inset of Fig. [TOjon a doubly- logarithmic scale clearly in- 
dicate that the algebraic 1/r behavior is observed, if at 
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Figure 9: Top: Displacement components Ux and Uy of 
the Iz interstitial as a function of distance r along the 
X-axis and y-axis, respectively (solid lines). Also plotted 
are the displacements computed from continuum theory 
by Ewald summation according to Equ. (|46)) (dashed 
line), and simple 1/r-behavior (dotted line). The inset 
shows the region close to the defect location. A defect 
strength of 7 = 0.2347ct^ was used. Bottom: Displace- 
ment components as above for an Id interstitial with 
7 = 0.2425cr2. 



all, only for large system sizes and in a limited distance 
range. 

A comparison of the Ewald summation results with 
the numerical calculations over the whole simulation cell 
is shown in Fig. 1111 The color coded map represents the 
relative deviation (see Equ. of Up(r) from Uc(r). In 

the figure, regions of large and small relative deviation 
are colored in red and blue, respectively. We find that 
the I2 and 1^ configuration show a relative deviation 
between 1% to 5% over the whole range. Only in the 
core region of the defect the deviations are larger. For 
the Id defect the deviations are larger and between 10% 
to 20% also far away from the defect. These deviations 
are due to discrepancies both in orientations as well as 
magnitude. 

The energy density of a point defect with periodic 
boundaries, calculated from Equ. for the displace- 
ment field of Equ. is depicted in Fig. [T^l As for 
the point defect in a circular rigid container discussed 
in Sec. 13.3.11 the energy density becomes constant for 
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Figure 10: Displacement component function 
of distance r from the defect along the x-axis obtained 
from simulations (solid lines) and according to Equ. (j46p 
(dashed hues) for system sizes = 26 x 30, 52 x 60, 78 x 
90, 104 X 120, 130 X 150, 156 x 180, 182 x 210, 208 x 240. 
Also shown is the 7/r behavior expected in an infinitely 
extended material (dotted line). The same value of 7 = 
0.23165(7^ was used in all cases. The vertical dotted 
lines indicate the distances of the cell boundaries from 
the origin for the various system sizes. In the inset the 
same curves are displayed on a logarithmic scale. 



large distances. This constant term arises from the work 
performed by the defect against the external pressure. 

4.2 Vacancies 

In the system studied in this paper, vacancies gener- 
ate displacement patterns that are considerably more 
intricate than those of interstitials, as can be inferred 
from a comparison of Figs. [T] and [21 While in the 
case of interstitials the displacement vectors essentially 
point away from the defect site, vacancies have displace- 
ment fields which point outward or inward depending 
on the position relative to the defect. For instance, in 
the V2 vacancy shown in Fig. [2^ the displacement vec- 
tors point towards the defect site along the x-axis, but 
away from the defect along the y-axis. Between the two 
axes, vortex like structures occur. A similar displace- 
ment pattern with alternating displacement directions 
forms also for the Va vacancy shown in Fig. [2j:. This 
behavior observed in the core region around the defect 
can not be reproduced by the simple defect model used 
here for the continuum theory calculations. The dis- 
placement field obtained in this model is either oriented 
towards the defect or away from it depending on the 
sign of the defect strength 7. For the V3 vacancy, on the 
other hand, all displacement vectors point inward and 
no such complications occur. Nevertheless, the displace- 
ment magnitudes shown in Fig. [T2]for the three vacancy 
configurations follow qualitatively the form predicted by 
continuum theory. For all three configurations, varying 
the defect strength 7 can lead to better agreement in 
particular directions (e.g., along the x- or y-axis), but 
not over the entire plane. 
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In Fig. [14] we depict the relative deviations of the 
displacement of the particle simulation from the contin- 
uum theory calculated according to Equ. [43l Even far 
from the defect the displacement field obtained from the 
particle simulation does not become isotropic such that 
it cannot be reproduced by the continuum theory cal- 
culations. The best agreement is found for V3 (Fig. [Ml 
center), in which case the anisotropy of the displacement 
field is less pronounced. The failure of the point defect 
model to reproduce the displacement patterns of vacan- 
cies, however, does not imply that elasticity theory is 
unsuitable for the description of such defects. Rather, 
the point defect model used here, which consists of two 
orthogonal pairs of opposing forces, appears to be to 
simple to capture the complex displacement patterns in- 
duced by vacancies. 



homogeneous neutralizing background that is seamlessly 
incorporated in the Ewald sum solution. The neutral- 
izing background satisfies the condition of fixed volume 
by exactly compensating for the volume change caused 
by the introduction of the point defect. The volume 
compensation leads to an additional term in the energy 
density related to the work done by the defect against 
the external pressure. Depending on the pressure, this 
energy can contribute significantly to the total defect 
energy. 

While for interstitials the elasticity theory calcula- 
tions carried out for a simple point defect model lead to 
good agreement with the particle calculations in the core 
region around the defect, large deviations are observed 
for vacancies. These discrepancies are due to the more 
complex displacement patterns of vacancies and better 
defect models are required to capture this behavior. 



5 Conclusion 

Interstitials and vacancies occur in various configura- 
tions generating displacement fields with symmetries 
that differ from the symmetry of the underlying triangu- 
lar lattice. Near the defect, the displacement fields are 
highly anisotropic and strongly dependent on the atom- 
istic details of the interactions. In this distance regime, 
linear elasticity theory brakes down due to discrete lat- 
tice effects and non-linearities of the potential. For dis- 
tances larger than about 10-15 lattice spacings, however, 
elasticity theory is valid. To establish this validity, it is 
crucial that corresponding boundary conditions are used 
both in the continuum calculations and the particle- 
based numerical simulations. Since simulations are usu- 
ally carried out with periodic boundary conditions in 
order to minimize finite size effects, the same boundary 
conditions must be employed also in the continuum cal- 
culation. If different boundary conditions are used, the 
long-range nature of elastic displacement fields can lead 
to considerable discrepancies even at length scale where 
elasticity theory is expected to hold. 

In this paper we have formulated the elastic the- 
ory problem in a way that makes it formally identical 
to the problem of determining the potential of a point 
charge in electrostatics. While here we have focused 
on two-dimensional systems, the same formalism ap- 
plies also to three dimensions. Under periodic boundary 
conditions, this two-dimensional electrostatics problem 
has been solved using the method of Ewald summation 
[m [42] , in which the solution is expressed in terms of 
two rapidly convergent sums, one in real space and one 
in reciprocal space. 

The solution of the electrostatics problem can be sim- 
ply transferred to the continuum theory of the point 
defect. In this case, the role of the charge density in 
electrostatics is played by the dilatation, i.e. the local 
relative volume change. Accordingly, the charge neu- 
trality required by the periodic boundary conditions in 
electrostatics corresponds to the condition of fixed vol- 
ume in the elasticity theory. This requirement leads to a 
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Figure 12: Energy density e of an I2 interstitial in a pe- 
riodic box as a function of the distance r measured along 
tlic x-dircction (solid line) and the y-direction (dotted 
line). Also shown is the energy density obtained from 
continuum theory (dashed line) according to Equ. 
for the displacement of Equ. (|46)) . 



Figure 11: Color coded relative deviations of the dis- 
placement Up(r) obtained numerically from the contin- 
uum theory prediction Uc(r) as a function of x and 
y calculated according to Equ. (j43|l . From top to 
bottom, the figures depict the relative deviations for 
the I2 interstitial (7 = 0.23165cr^), the /a interstitial 
(7 = 0.2347(7^), and the h interstitial (7 = 0.24250-^). 
The whole simulation cell of dimensions L,j. — 215.28(t 
and Ly = 215.12cr is shown. Colors are assigned on 
a logarithmic scale which runs from 10""^ (blue) to 10^ 
(red). The white, black, and blue contour lines represent 
a relative error of 1%, 2% and 5%, respectively. 
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Figure 13: From top to bottom: Absolute value of the 
displacement components (black solid line) and Uy 
(black dashed line) of a V2, V3 and Va vacancy. Also 
plotted are the absolute values of the predicted displace- 
ment field from the Ewald summation Equ. ([46]) for the 
defect strengths 7 = -I.TSct^ {V2), 7 = -0.29ia^ {V3), 
and 7 = — 1.89cr^ i^a), found by minimizing the relative 
error of Equ. l43l over the whole xy-plane (dashed line), 
as well as l/r behavior (dotted line) . 




Figure 14: Color coded relative deviations Equ. (|43p of 
the displacement obtained numerically for three different 
vacancy configurations from the prediction of continuum 
theory as a function of x and y. Color code and system 
size is the same as in Fig. ([TT]) . From top to bottom: 
V2 vacancy (7 = -l.TScr^), V3 vacancy (7 = -0.294ct2), 
and Va vacancy (7 = — 1.89cr^). 



16 



